Limits of quotients of real analytic functions in 

two variables 



Carlos A. Cadavid*, Sergio Molina**, Juan D. Velez**, 

* Corresponding Author, Universidad EAFIT, 
Departamento de Ciencias Basicas, 
Bloque 38, Office 417. Carrera 49 No. 7 Sur -50, 
Medellm, Colombia, ccadavid@eafit.edu.co, (57) (4)-2619500. Fax (57) (4)-3120649. 
** Universidad Nacional de Colombia, Departamento de Matematicas, 
Calle 59 A No 63 - 20 , Oficina 43-106, Medellm, Colombia, 
sdmolina@gmail . com , j dvelez@unal . edu . co 

November 9, 2010 



Abstract 

Necessary and sufficient conditions for the existence of limits of 
the form 

hm — 

(x,y)->(a,b) g{x,y) 

are given, under the hipothesis that / and g are real analytic functions 
near the point (a, b), and g has an isolated zero at (a, b). An algorithm 
(implemented in MAPLE 12) is also provided. This algorithm deter- 
mines the existence of the limit, and computes it in case it exists. It is 
shown to be more powerful than the one found in the latest versions 
of MAPLE. The main tools used throughout are Hensel's Lemma and 
the theory of Puiseux series. 

Keywords: Limit, Real Analytic Function, Hensel's Lemma, Puiseaux 
Series 
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1 Introduction 



In the usual calculus courses one is asked to determine the existence of limits 
of the form 

hm — 

(x,j/)-Ka,&) g{x,y) 

where / and g are real analytic functions (typically, polynomials or trigono- 
metric and exponential functions) defined in an open disk centered at a point 
(a, b) in M 2 . The standard strategy for solving this problem consists in study- 
ing the existence of the limit along various simple trajectories, such as straight 
lines, quadrics, cubics, etc., with the hope that either one of them fails to 
exist or two of them differ. If they all coincide, then one tries some other 
ad hoc trajectories. If all that fails, one tries to prove its existence by some 
theoretical methods. 

In this paper we develop a theoretical method which completely solves this 
problem. An algorithm for polynomials based on this method is implemented, 
which proves to be more powerful than other existing routines. 

An application of Weierstrass' Preparation Theorem allows to reduce the 
problem to the case where / and g are monic polynomial functions in the 
variable y, whose coefficients are real series in the variable x. Next, a discrim- 
inant real curve is constructed using Lagrange Multipliers with the property 
that the limit exists, if and only if, it exists along this curve. Then Hensel's 
lemma, some Galois theory, and the theory of Puiseaux series are used to 
parametrize the various branches of the discriminant curve and select the 
real ones. All these steps are done in a constructive manner making it pos- 
sible to implement this method in an algorithmic way. In this article an 
algorithm was implemented for polynomial functions. 

2 Theory 

2.1 Reduction to the case where / and g are polyno- 
mials 

After a translation, we may assume that (a, b) is the origin. 

Let us denote by S = R {x, y} the ring of power series in the variables 
x, y with real coefficients having positive radius of convergence around the 
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origin. If h(x, y) belongs to S, the order of h in the variable y is defined to 
be the smallest integer r such that h(y) = h(0, y) has the form 

h(y) = a r y r + a r+ iy r+ + • • ■ , with a r ^ 0. 

(If h(y) — the order is defined to be +00.) 

It is not difficult to show that given hi, . . . , h n in S — {0} there exists 
an integer v > 1 such that after a change of coordinates of the form x' = 
x + y v ,y' = y, each series h'^x' , y') = hi(x + y v ,y) is of finite order in the 
variable y' |GLS] . The essential tool for the reduction is the following lemma. 

Lemma 1 (Weierstrass) Let h be an element of S = R {x, y} of order 
d in y. Then there exists a unique unit u(x,y) G S and unique real se- 
ries ai(x) , . . . , ad(x) with positive radii of convergence such that h(x,y) = 
u(x, y)(y d + ai{x)y d - 1 + ■■■ + a d {x)) iGLSf . 

Since the existence of the limit and its value is obviously independent 
of the particular choice of local coordinates, we may assume that f(x,y) = 
u(x,y)fi(x,y) and g(x,y) = v(x,y)gi(x,y), where u{x,y) and v(x,y) are 
units and 

fi = y d + ai{x)y d ~ l + ■ ■ ■ + a d {x) 
9i = V b + ci{x)y b - 1 + • • • + c b {x) 

are monic polynomials in R {x} [y] . Since units do not affect the existence 
of the limit, there is no loss of generality in assuming also that u and v are 
equal to 1. 



2.2 Discriminant variety for the limit 

The following proposition provides a necessary and sufficient condition for 
the existence of the limit. 

Let / = y d + ai(x)y d ~ l + ■ ■ • + a^(x) and g = y b + ci(x)y b ~ 1 + ■ ■ ■ + Cb(x) be 
monic polynomials in R{x} [y] , and D = D p (0) C R 2 a closed disk centered 
at the origin with radius p > 0, such that each ai(x), Cj(x) is convergent in 
D. Let us denote by h' the polynomial ydq/dx — xdq/dy, where q denotes 
the quotient q = f ' / 1 g, and by h" (the numerator of h') the polynomial 

h" = y (g-L - f?l) - X (a— - f— \ (1) 

\ dx dx J \ dy dy J 
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Let X be the variety cut by h' in the puncture disk, i.e., 



X = {(x, y) G D : (x, y) ^ (0, 0) and h"(x, y) = 0} . 
With this notation we have the following proposition. 
Proposition 2 Let q(x,y) = f(x,y)/g(x,y). The limit 

lim ?/) 

exists and equals L £ R, if and only if for every e > there is < 5 < p 
such that for every (x,y) G X fl D$, the inequality \q(x,y) — L\ < e holds. 

Proof. The method of Lagrange multipliers applied to the function q(x,y), 
subject to the condition x 2 + y 2 = r 2 where < r < p, says that the extreme 
values taken by q(x, y) on each circle C r (0), centered at the origin and having 
radius r, occur among those points (x,y) of C r (0) for which the vectors 
(dq/dx, dq/dy) and (x, y) are parallel, which amounts to ydq/dx — xdq/dy = 
0. Let us assume that given e > there exists < 5 < p such that for every 
(x,y) G X fl D$ the inequality \q(x,y) — L\ < e holds. Let (x, y) G D$ and 
r = a/x 2 + y 2 . If ti(r),t 2 {r) G C r (0) are such that q(h(r)) = min teC v(o) q{t), 
and q(t 2 (r)) = max (eCr(0) q(t), then 

q(tx(r)) - L < q(x, y) - L < q(t 2 (r)) - L, 

for every (x,y) G C r (0). Since ti(r),t 2 (r) G X fl D$ we have 

-e < q{h(r)) - L and q(t 2 (r)) - L < e 

and therefore \q(x, y) — L\ < e, for all (x, y) G D$. 

The "only if part is immediate from the definition of limit. ■ 

2.3 Hensel's Lemma 

Let us fix an integer n > 0. A linear change of coordinates of the form 
fi(x, y) = f(x+ny, —nx+y), and gi(x, y) = g(x + ny, —nx + y) does not alter 
the limit of the quotient as (x, y) approaches the origin. It is easy to see that 
this change of coordinate transforms (JT|) into a monic polynomial multiplied 
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by a nonzero constant. By the chain rule we have that h"(x + ny, —nx + y) 
is equal to 

(-nx + y)[g 1 {df/dx) 1 - f^dg/dx)^ - (x + ny)[g l (df/dy)i - f^dg/dy)^, 

where (df/dx)i(x,y) = df/dx(x+ny,—nx+y), and similarly with (df/dy)i, 
(dg/dx)i, (dg/dy)i. We will denote h"(x + ny, —nx + y) simply by h(x,y). 

Our next goal is to parameterize the curve h(x, y) = 0. For this purpose 
we will use Hensel's Lemma ([E]). Let us denote by k a arbitrary field, by 
R the ring of formal power series in the variable x, with coefficients in k, 
R = k[[x]}, and by k((x)) its field of fractions. R is a local ring (R, m) whose 
maximal ideal is m = (x). For each h(x,y) e R[y] monic in the variable y, 
let us denote by h its reduction modulo m, i.e., h = h(0,y). 

Lemma 3 (Hensel's Lemma) Let F{x,y) be an element of R[y] monic in 
y, and let us assume that F = gh is a factorization in k[y] whose factors are 
relatively prime, and of degrees r and s. Then there exist unique G and H in 
R[y] with degrees r and s, respectively, such that: 

1. G = g and H = h 

2. F = GH 

In order to construct a parameterization of h(x,y) = 0, Puiseaux series 
are used which we review next. 

2.4 Puiseaux Series 

Let us denote by L the quotient field of fractions of R = C[[x]}, which consists 
of Laurent series. Let L be an algebraic closure of L. For each positive integer 
n we will denote by x l ^ n a fixed n-th root of x in L. It is clear that the n-th 
roots of x are 

6x 1/n ,e 2 x 1/n ,... 1 n - 1 x lln ,x lln 

where 9 is any primitive n-th root of unity. It is easy to see that the poly- 
nomial t n — x is irreducible in L[t] and therefore L C L(x 1 ^ n ) is an extension 
of degree n. Consider the directed system consisting of the positive natural 
numbers (partially) ordered by divisibility, i.e., n < m if and only if n\m. 

The direct limit lim ^ L(x l l n ) will be denoted by L*. This limit can be 

identified with the field U n L(x 1//n ) C L. Each element a of L* can therefore 
be written in the form a = CkX qk , with Cf. G C, and exponents G Q such 
that: 
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1. q 1 < ■■■ <q r < 



2. There is an integer b so that each exponent can be written as qi = a,i/b, 
for some integer a,. 

The least exponent in the expression for a, qi, is called the order of a. It 
is a well known theorem that given a monic polynomial 

h(x, y) = y d + h^y^ 1 + • • • + h d (x) 

in R[y], there is an integer N > such that h can be factored completely in 

L(x l l N ) C L* as 

h(x, y) = (y- a^x 1 '")) ■■■{y- a^x 1 ^)), (1) 

where each <7j(t) is an element of C[[i]], i.e. a formal power series. Moreover, 
it can be seen that this series has positive radius of convergence and therefore 
defines a holomorphic function (cf. [GLS] ). Using this result, it is possible 
to parameterize the curve 

X = {{x, y) G C 2 : {x, y) ^ (0, 0) and h{x, y) = 0}. 

The proof of the existence of ([1]) can be done constructively using Hensel's 
Lemma (Lemma [H]) making it possible to determine which one of the series 
CTj(t) has only real coefficients. Such series will be called throughout, a real 
series. This in turn allows us to parameterize each one of the trajectories in 
Xnl 2 which go through the origin. It will be shown that these are the only 
ones that are relevant, since for a holomorphic function to be real valued on 
a real sequence approaching zero, it must have a series expansion around the 
origin with only real coefficient. 

The parameterization of the zeroes of h can be done by observing that 

h( xN ,y) = (y - •■■(y- o^O)) 

and consequently X = {(x,y) G C 2 : h(x,y) = 0} is the union of the sets 
Xi = {(z , &i(z)) : z G C}. This allows us to prove the following central 
result. 

Theorem 4 Let ci(z), . . . , &i(z), I < d, be the real series in the equation (Q]] 
which go through de origin (i.e. <7j(0) = 0, for i = 1, . . . , /) . Then the limit 

hm — r 

(x, y )->(p,o) gi(x, y) 
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exists if and only if 

«-» 9i (t N , o,(t)) 
exists, for i — 1, . . . , I, and L\ = ■ ■ ■ = L\. 

2.5 Newton's automorphism 

For each rational number g ^ there exists a homomorphism a q : L* — > 
L*, which sends x to x q and fixes the subfield C. This homomorphism is 
constructed by first defining a homomorphism from C[x] into L* which sends 
x to x Q , then extending it to C[[sc]], and then to the field of fractions C ((#)). 
Since L* is an algebraic extension of L, this homomorphism extends to a 
homomorphism a q from L* to L. It is clear that the image a q lies inside 
L*, and therefore one can regard a q as an endomorphism of L*. For p ^ 
rational, let (3 qiP : L* [y] — > L* [y] be the extension of a 9 obtained sending y to 
j/x p . It is clear that (3 gtP is invertible and its inverse is /3i/ q - p / q . 

With these preliminaries we can now state the following fundamental 
theorem ([M])- Even though this result is well known in the literature, we 
provide a "constructive" proof, since it is the very heart of the procedure sus 
in the algorithm limite, whose code we give at the end. 

Theorem 5 Every polynomial h = y d + h 1 (x)y d ~ 1 + • • • + h^(x) with coeffi- 
cients hi(x) in L* can be factored into linear factors h = (y — ax) • • • (y — a d ), 
with (jj G L* . Even more, if each hi(x) belongs to C[[x]], then there exists a 
positive integer n such that all o~i E C[[x 1//n ]]. 

Proof. The proof proceeds by induction on the degree of h. The homo- 
morphism <f) : L*[y] — > L*[y] sending y to y — (hi(x)/d), and fixing every 
element of L*, is invertible and its inverse is the homomorphism that fixes 
each element of L* and sends y to y+ (hi(x)/d). A simple calculation shows 
that <j)(h) is a polynomial with coefficients in L* such that the coefficient of 
y d ~ l is zero. Then 

<P(h)=y d + b 2 (x)y d - 2 + --- + b d (x). 

Let us denote by U{ the order of hi (x) , and let 2 < r < d be the least index for 
which u r /r = mm{ui/i : 2 < % < d}. Let us define if) = fi r ,u r '■ L*[y] — > L*[y}. 
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ip sends x to x r , y to x Ur y (and its inverse sends x to x 1 ^, and y to x Ur ^y). 
Clearly, 

V#(/i)) = x dUr y d + b 2 (x r )x (d - 2)Ur y d - 2 + ■■■ + b d (x r ) (2) 

= x dUr {y d + x- 2ur b 2 (x r )y d - 2 + ■■■ + x' kUr b k (x r )y d - k + ■■■ + x~ dUr b d {x r )). 

(3) 

The order of each term x~ kUr bk{x r ) of the polynomial inside the parentheses 
is given by — ku r + ruu > (since ut/k > u r /r). Furthermore, the order of 
x~ rUr b r (x r ) is 0. Therefore, if 

F' = y d + x- 2Ur b 2 (x r )y d ~ 2 + ■■■ + x~ kUr b k {x r )y d - k + ■■■ + x- dUr b d (x r ), 

by taking N large enough, we have that F = F'(x N ,y) G C[[x] ][?/], and F 
admits a modulo x reduction / = F G C[y] having at least two distinct roots. 
For if / had a single root c, then it is impossible that c = because the order 
of x~ rUr b r (x r ) is zero. And if c ^ then / = (y — c) d = y d — dcy^ 1 + • • • 
would have a nontrivial term in which is also impossible. Therefore, 

since / G C[y], one obtains / = /i/ 2 , with f 1 and / 2 monic and of degrees 
strictly smaller than d. Hensel's Lemma guarantees the existence of a lifting 
F = FiF 2 with Fi and F 2 monic, and of degrees strictly smaller than the 
degree of F. Consequently, F'(x,y) = F[F 2 with F( = Fj(x 1//7V , y). In 
conclusion, ip((j)(h)) = x dUr F[F 2 . By the induction hypothesis we know that 
F{Fj = rij=i(y ~ °j)> °j e ^* an< i therefore 

MO = a**t r il>-\F[F*) (4) 

= n(y - x^'r^i)) (6) 

is also a product of linear factors. This finishes the induction. The last claim 
in the theorem also follows by induction on the degree of h. In fact, if the 
coefficients of h belong to C[[x]], then F' G C [[#]][?/] and consequently F[ and 
F 2 also belong to C[[x]} [y]. The induction hypothesis guarantees the existence 
of a positive integer m such that F[F 2 = Yl^ =1 (y — (Tj), with <jj G C[[x 1,/m ]]. 
Furthermore, in the polynomial (j>{h) all the u k are nonnegative and therefore 
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u r /r > 0. Consequently, each element x Ur ^ r ip — x Ur l r o-j(x 1 l' r ) belongs 

to qix 1 /"]] with n = rm, and therefore the desired factorization for h is 
obtained. ■ 

It easily follows from the last part of the proof of theorem [5] shows that 

Remark 6 If h = y d + hi(x)y d ~ 1 + • ■ ■ + hd(x) is a monic polynomial with 
coefficients mC[[x]], then there exists a power r > 0, and polynomials gi(x,y) 
and g 2 (x,y) in C[[x]][y], which are monic in the variable y and of degrees 
di, d 2 < d, such that h{x r , y) = g±(x, y)g 2 (x, y). 

Theorem |5] admits the following refinement (JM]). 

Theorem 7 Let h = y d + /^(x)?/^ 1 + • • • + hd{x) be a monic polynomial 
with coefficients in C[[x]] that is irreducible over L = C((x))[y]. Then, if uo 
denotes a primitive d root of unity, there exists a(t) = Yl'kLo ^ such that 



The following result allows to identify the real series in the factorization 
given in ([5]). 

Lemma 8 Let F = y d + bi(x)y d ~ 1 + ■ ■ • + b d (x) be a monic polynomial of 
degree d in the variable y and whose coefficients are real power series, i.e. 
bi(x) E JR[[x]] . Then F = (y — r) d , with r 6 M ; if and only if its factorization 
in Puiseaux series has the form 



h = (y — a(ujx 



!/«*)) .-.(y- a(iJ d X 1/d )) 



where a(u r x 



/d ) = Y.k=ocMx 1 ' d ) k . 



s 



F^F t 



(7) 



i=l 



where 
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Proof. Let us first consider the only if part of the equivalence. We proceed 
by induction on d. If the degree of F is d = 1, then F = y — b\(x), with 
b\(x) G IR[[x]], is a factorization in Puiseaux series and 61 (0) = r. If d > 1, by 
remark [6], there exists an integer N > such that F(x N , y) = G(x, y)H(x, y) 
where G and H are polynomials which are monic and of degrees d\ , di < d 
in the variable y. Thus 

F = F(0,y) = G(0,y)H(0,y). 

And this implies that 

G = (y - r) d \ and H = (y - r) d \ 

By Hensel's Lemma applied to the ring and in particular, by its 

claim about uniqueness, one obtains that G and H are polynomials with 
coefficients in By the induction hypothesis G and H can be factored 

in the form ([7]), so F can also be factored in this way. Conversely, if ([7j) 
holds, then it follows, by taking iV = d\ . . . d s , that 

s 

F(* N ,y) = 11^ - °& ei ))(v ~ °i{ui* ei )) ■ ■ ■ (y - ^t 1 ^)) (8) 

where = N/di and cjj(t) is a real series. By replacing x by in (JSJ) one 
obtains 

F = F(0, y) = Y[(y - ^(0))* = (y - r) d . 

i=l 

■ 

Let now F = y d + b\(x)y d ~ 1 + • ■ ■ + bd(x) be a monic polynomial whose 
coefficients are in K[[x]] and let us denote its reduction modulo x by /. We 
can write 

f = (y- nfi ■■■{y- r s ) d 's(y - Cl ) d ^y - cl) d ^ ■ ■ ■ (y - Cl )\y - ci) d \ 

where r 1; . . . ,r s are the real roots of /, and Cj,Q are the nonreal ones. Let 
us define fi(y) — (y — ri) di and 

9i {y) = {y- Cl t{y - = (y 2 - a t y + ^) d % 

with «j, fa real. Hensel's Lemma provides us with a lifting of the factorization 
fx--f»9\-'9h of the form 

F — Fx ■ ■ ■ F s Gi ■ ■ ■ Gi 
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i.e., Fi = fi and G, = g^. From the proof of Hensel's Lemma it follows that 
each Fi is a monic polynomial in the variable y with coefficients in M[[x]]. 
Each Gi admits a factorization Yif=i irreducible factors in C[[a;]][y], 

and, by Gauss' Lemma, also in C((x))[y] (see [L]). Notice that if Gij has 
degree e^ , then Ylf=i e «i = 2c?i- By theorem [7] each can be factored as 

Gij 

Gi^lKy-a^x 1 ^)), 

k=l 

where Uj is an efj primitive root of unity and <Tj(t) £ C[[i]]. If we let = 
e»i . . . e iqi , it is clear that 

i=i fe=i 

and as a consequence 

^(o, 2 /)=n(y- ( T,(o)) e -. 

It follows that qi = 2 and cxi(O) = q and a"2(0) = Cj (or the other way round). 
Therefore none of the <jj is a real series. By lemma El F k = YYh=i Fkh, and 

F kh = {y - a hk ^ 1/dhk )){y ~ <r hk {u hk x l ' d ^)) ■ ■ ■ (y - a^f^x 1 '^)), (9) 

with <Jhk{t) being a real series. The following theorem summarizes what has 
been achieved so far. 

Theorem 9 Let F = y d + bi(x)y + • • • + bd(x) be a polynomial that is 
monic in the variable y and whose coefficients lie in and let f be its 

reduction modulo x. Then f can be written as 

/ = (y - n)* ■■■{y- r s ) d '^y - Cl ) d ^y - cl) dl • • • (y - d)*(y - ci) d K 

Let 

fi{y) = (y~ Ti) d K 9i(y) = (y- c)*(j/ - i%* = (y 2 - a t y + ft)* 

with ri, «j, [3i real. Hensel's Lemma gives a lifting of the factorization of f = 
fi ■ ■ • fs9i ■•■9i 

F — Fi ■ ■ ■ F S G\ ■ ■ ■ Gi 

with Fi = fi and Gi = gi. Then, in the Puiseaux series factorization of F 
the only real series occur in the decomposition into linear factors of the Fi. 
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Thus if X denotes the curve in C 2 formed by the zeroes of h then 
Ifll 2 = U* =1 {(t d , ai(t mi )) :te[a,b]G R}. 

3 Algorithm for the computation of limits 

As an application of the results established in the previous section, we present 
in this section an algorithm implemented in Maple 12 for the determination 
of limits of the form 

hm — 

(x,y)->(o,o)g{x,y) 

where /, g e M[x, y] and g has an isolated zero at (0, 0). The complete routine 
is called limite and comprises several subroutines that are documented next. 

The routine repeticion takes as input a list L and forms a list of lists, 
each one formed by each element of L, repeated as many times as it appears 
in L. For instance, if L — [1, 2, 1, 3, 1, 2, 4, 3] then repeticion produces the list 
[[1, 1, 1], [2, 2], [3, 3], [4]]. The Maple 12 code for this routine is as follows. 

> repeticion: =proc(L) 

> local S,j,H,i; 

> H : = convert ( convert ( L , set ), list ) ; 

> for i from 1 to nops(H) do S[i]: = []; 

> for j from 1 to nops(L)do 

> if L[j]=H[i] then S[i]:=[op(S[i]),H[i]]; end if; 

> end do; end do; 

> RETURN ( [seq(S [i] , i=l..nops(H))]); 

> end proc: 

The routine suprime takes as input a list L and eliminates its redundancy. 
For instance, if L = [1, a, 1, 3, a, b, c], then suprime returns G = [1, a, 3, b, c]. 
The code for this routine is 

> suprime: =proc(L) 

> local S,i,G; G:=[L[1]]; 

> for i from 2 to nops(L) do 

> S:=convert(evalf(G),set); 

> if (evalb('in'(evalf(L[i]),S)) = false) then G:=[op(G),L[i]]; 

> end if; 

> end do; 
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> RETURN(G); 

> end proc: 

The routine poli takes as input a list L whose elements are complex 
numbers and constructs another list containing a polynomial of the form 
(y — r) d for each real r that appears exactly d times in L, and a polynomial 
of the form {{y — z)(y — z)) d for each nonreal z appearing together with its 
conjugate z exactly d times in L. For example, if L — [1, 2 — i, 1, 2 + i] then 
poli returns [(y — l) 2 , (y — (2 — i))(y — (2 + i))\. Its code is 

> poli:=proc(L) 

> local g,H,i; 

> H:=repeticion(L); 

> for i from 1 to nops(H) do 

> if Im(evalf(H[i][l]))=0 then g[i]:=(y-H[i][l])^nops(H[i]); else 

> g[i]:=((y-H[i][l])*(y-conjugate(H[i][l])))^nops(H[i]); end if; end do; 

> supr ime ( [seq (g [i] , i= 1 . . nops ( H ) ) ] ) ; 

> end proc: 

The routine mochar takes a polynomial f(x, y) = a (x)y d +- ■ - + ak(x)y k + 
■ ' ' + a d{ x ) and eliminates from each coefficient those powers of x which are 
larger than n, i.e. it calculates f(x, y) modulo x n+1 . The code for this routine 
is 

> mochar: =proc(f,n) 

> local c,d,i,g; 

> d:=degree(f,y); g:=0; 

> for i from to d do 

> c[d-i]:=mtaylor(coeff(f,y,d-i), [x], n+1); 

> g:=g+c[d-i]*y"(d-i); 

> end do; 

> collect (g,y); 

> end proc: 

The routine monico takes a polynomial / in the variable y and divides it 
by the coefficient of the highest power of y. Its code is 

> monico: =proc(f) 

> local c,d; 

> d:=degree(f,y); 

> c:=coeff(f,y,d); 
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> collect(expand(l/c*f),y); 

> end proc: 

The routine Hensel has four entries. The first entry is a polynomial 
F(x,y) en C[x][y], which is monic in y. The second and third entries are 
polynomials g(x), h(y) such that F(0,y) = g(y)h(y). The fourth entry is 
an integer n. Hensel calculates polynomials G(x, y) and H(x, y) such that 
F = GH modulo x n+1 . The code for this routine is 

> Hensel:=proc(poly,gg,hh,n) 

> local L,H,G,i,l,f,g,h,t; 

> f[0]:=coeff(poly,x,0); 

> g[0]:=gg; 

> h[0]:=hh; 

> f[l]:=coeff(poly,x,l); 

> gcdex(g[0],h[0],f[l],y,'s','t'); h[l]:=s;g[l]:=t; 

> for i from 2 to n do 

> f[i]:=coeff(polyx,i);l[i]:=f[i]-sum(g[j]*h[i-j],'j'=l..i-l); 

>gcdex(g[0],h[0],l[i],y,'s','t'); 

> h[i]:=s;g[i]:=t; 

> end do; 

>H:=sum(h[j]*x"(j),'j'=0..n); 

> G:=sum(g[j]*x"(j), , j'=0..n); 

> L: = [mochar(G,n), mochar(H,n)]; RETURN (L[1],L [2]); 

> end proc: 

The routine henselgen takes as entry a polynomial f(x, y) which is monic 
in y, a list L of polynomials fi(y), all monic in y, pairwise relatively prime and 
such that /(0, y) = fi(y) ■ ■ ■ f r (y), and an integer n > 0. This procedure re- 
turns polynomials G±(x, y), G r (x, y) such that f(x, y) = Gi(x, y) • • • G r (x, y) 
modulo x n+1 , and Gi(0,y) = fi(y). The Maple 12 code for this routine is 

> henselgen: =proc(f,L,n) 

> local P,G,H)T,S,i; S:=convert(L,set); H:=[]; 

> for i from 1 to nops(L) do 

> T:= convert(S minus {L[i]},list); 

> P [i] : =expand(product (T [k] ,k= 1 . .nops(T) ) ) ; 

> G[i]:=Hensel(expand(f),expand(L[i]),P[i],n)[l]; H:=[op(H),G[i]]; 

> end do; 

> RETURN (H); end proc: 
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The routine orden receives a polynomial "poly" 

f(x,y) = y d + b 1 (x)y d - 1 + --- + b d (x), 

with each bi(x) being a Laurent polynomial in x. It computes u = min{u[i] : 
i = 1, ...,d}, with u[i] being the order of bi(x) (i.e. the degree of the least 
degree term in and returns [r, u[r]], with r the smallest index i such 

that u[r]/r = mm{u[i]/i}. Its code is 

> orden: =proc(poly) 

> local u,i,b,r,m,d; 

> d:=degree(poly, y); 

> b[l]:=coeff(poly,y,d-l); u[l]:=ldegree(b[l]); r:=l; m:=u[l]/r; 

> for i from 2 to d do 

> b[i]:=coeff(poly,y,d-i); u[i]:=ldegree(b[i]); 

> if (u[i]/i < m ) then 

> m:= u[i]/i; r:=i; 

> end if; 

> end do; RETURN([r,u[r]]); 

> end proc: 

The routine sus receives a polynomial "poly", 

poly(x, y) = y d + h^y*" 1 + ■■■ + b d (x), 

whose coefficients are Laurent series in x and computes: 

1. A polynomial g (x, y) in C((x))*[y] satisfying tp(f)(f(x, y)) = x du ^g(x,y), 
(for notation see |5]) where / denotes the polynomial that is obtained 
from poly after performing the substitution y — y — b x (x) /d, i.e. 

f(x,y) =poly(x,y-b 1 (x)/d) = y d + c 2 (x)y d ~ 2 H bc d (x). 

Then r is equal to the value of the index i where the minimum value 
of : u[i] = degree of Cj(x)}, is attained, and u[r] is equal to the 

order of c r (x). where 

^ : C((x))*[y] ^ C((x))*[y], 

denotes the automorphism obtained by composing the map <fi : f(x,y) i— 
f(x,y- 6i(x)/d), with if) : f(x,y) i — > f(x r } yx u[r] ). 
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2. sus returns the triple [g,r, u[r]], with g(x,y) = x du ^ip(j)(f /x, y). 
Warning: 

1. If u[r] = oo ( which happens in case f(x,y) = y d , or equivalently if 
poly(x,y) — (y — bi(x)) d ), then sus returns the triple [(y — bi(x)) d , 1, oo]. 

2. Notice that r and u[r] correspond to the polynomial f(x,y) obtained 
from "poly" after performing the linear substitution that eliminates the term 
of degree d — 1 in y, but not to the polynomial "poly" itself. 

The code for the routine sus is 

> sus: =proc (poly) 

> local g,q,f,b,d,u,r; 

> d:=degree(poly,y); 

> b[l]:= coeff(poly, y d-1); 

> f:=collect(simplify(subs( y=y-b[l]/d, poly)),y); 

> r:=orden(f)[l]; u[r]:=orden(f)[2]; if u[r]= infinity then RETURN([factor(poly),r,u[r]]);end 

if; 

> q:=collect( expand(x" (-d*u[r])*subs( {x=x"r, y=y*x"u[r]}, f )), y); 
RETURN ( [q,r ,u [r] ] ) ; 

> end proc: 

The routine invsus takes as input an integer d, a polynomial b = b(x), 
integers r, u and a polynomial g(x, y) = y d + c\(x)y d ~ l + • • • + Cd(x). If ip -1 
is the automorphism determined by x i — > x 1 ^, y i — > yx~ u ^ r , and is the 
isomorphism defined as the linear substitution y i — > y + b(x) / d, then invsus 
computes x du ^ r (J)~ 1 ^~ 1 (g(x,y)). 

Observation: 

If d =grado(/(x, y)), b = b[l] the coefficient of y d ~ l in f(x,y), r the 
smallest index % such that u[r]/r = min{u[i]/i}, where each u[i\ is the degree 
of the coefficient of y d ~ l in the polynomial f(x, y — b/d) and u = u[r] the order 
of its r th coefficient c r (x), (quantities associated to f(x,y) in the procedure 
sus, and g(x,y) is a polynomial monic and of degree d in y, resulting from 
the previous procedure, i.e. such that g —sus(f(x,y))[l], and therefore 

g(x,y) = x- du[r] ip(j)(f/x,y)). 

Then invsus returns as a result f(x,y), because 

x du[r]/r ^0-V _1 (^" du[r V0(/ 7x,y)) = 

x du[r]/r x -du[r]/r ^-l^-l ^f/ x ^ y ^ = ffay). 
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Here is the code for sus. 

> invsus:=proc(d,b,r,u,g) 

> local D,t;D:=degree(g,y); 

> t:=root(x,r, symbolic); 

> simplify(t"(D*u)*subs({x=t, y=y*t"(-u)},g), symbolic); 

> collect (simplify (subs(y=y+b /d, %) ) ,y) ; 

> RETURN(%); 

> end proc: 

The routine reduction receives a polynomial f(x,y) which is monic in 
y, and an integer n > 0. If / is linear or has the form (y — bi(x)) d , 
the algorithm returns the triple [1, 1, f(x, y)]. Otherwise, it sets g(x,y) 
—sus(f) (hence ip<f>(f/x, y) = x du ^g(x, y)) and tries to verify whether g(0, y) 
has at least one real root. If this is not the case, the algorithm returns 
[1, 1, f(x, y)\. Otherwise, it factors g(0,y) in terms of the form g(0,y) = 
h\(y) ■ ■ -h m (y), where each hi(y) has the form (y — r,i) d % with r\ real or of 
the form hi(y) = [(y — z)(y — z)] d \ with z nonreal, lifts the factorization us- 
ing henselgen modulo x n+1 , and applies invsus in order to obtain an integer 
r > such that f(x r ,y) = qi(x, y) • • • q m (x, y). The procedure returns the 
list [r, [qi(x,y), ...,q m (x,y)]\. Here is the code for this routine. 

> reduccion:=proc(f,n) 

> local i,L,S,d,h,q,r,b,g,ra,H,u; 

> d:=degree(f,y); if d=l then RETURN([l,[l,f]]); end if; b[l]:=coeff(f,y,d- 
l);r:=sus(f)[2];u:=sus(f)[3];if u=infmity then RETURN([l,[l,factor(f)]]); end 

if; 

> g:=unapply(sus(f)[l],x,y); 

> if [fsolve(g(0,y),y)] = [] then RETURN([l,[l,f]]); end if; 

> S: = [solve(g(0,y),y)]; 

> L:=poli(S); 

> H:=henselgen(g(x,y),L,n); 

> for i from 1 to nops(H) do 

> h[i] :=unapply (invsus(d,b [1] ,r,u,H[i] ) ,x,y) ; 

> q[i] : =mochar (simplify (h[i] (x"r,y) , symbolic) ,n) ; 

> end do; 

> RETURN( [r, [seq(q[i] ,i=l . .nops(H) )]] ) ; 

> end proc: 

The routine fiscal checks the list L and returns the empty list [ ], if and 
only if all elements in L are powers of linear polynomials (i.e. each Li(x,y) 
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in the list L has the form (y — b 1 (x)) d ,d > 1) or each L^x, y) satisfies that 
sus(Li(x,y)) has only nonreal roots. Otherwise, it returns the same list L. 
Here is the code for this routine. 

> fiscal:=proc(L) 

> local ra,g,i; 

> for i from 1 to nops(L) do 

> g:=unapply(sus(L[i])[l],x,y); 

> ra:=[fsolve(g(0,y),y)]; 

> if (degree (L[i],y)>l and ra<>[] and sus(L[i]) [3] <> infinity ) then RE- 
TURN^); 

> end if; end do; RETURNQ]); 

> end proc: 

The routine cambio receives a list L = [r,[ri,r 2 , ...r i _x-, r i-, r i+\-, •••,' r n]], 
two integers i, m, and another integer b. The procedure returns a new list 
L = [br, [bri, br 2 , ...bri-i,ri, r,, r iy br i+ i, br n ] (from the i th position on it 
puts Ti repeated m times). The code for cambio is 

> cambio:=proc(L,i,b,m) 

> local n,S,h; n:=nops(L[2]); 

> if i=l then S:=[ L[l]*b, [seq(L[2][l] j=l..m),seq(b*L[2][j] j=i+l..n) ] ]; 
end if; 

> if i =n then S:=[ L[l]*b,[ seq(b*L[2][j],j=l..n-l),seq(L[2][n],j=l..m)] ]; 
end if; 

>S: = [L[l]*b,[seq(b*L[2][j],j = l..i-l),seq(L[2][i],j = l..m),seq(b*L[2][j],j=i+l. 

> S; end proc: 

The routine factoriza takes a polynomial f(x,y) that is monic in y and 
an integer n. The algorithm produces a list L = [LI, L2] with entries having 
the form LI = [fi(x, y), f n (x, y)\ and LI = [r, [ri, r n \\ such that r = 
T\ • • ■ r n and 

f(x r ,y) = f 1 {x r \y)---f n {x r - ) y), 

such that each fi(x,y) admits no real reduction, i.e. fi(x,y) is either of the 
form (y — bi(x)) d , with d > 1, or if gi(x,y) —sus(fi(x,y)), then gi(0,y) does 
not admit any real root. In other words, the only real Puiseux series in the 
factorization of f(x,y) are the ones given by fi(x T1 ^ r , y), f n (x Tn ^ r , y), and 
therefore if 

f(x, y) = (y-(7 1 (x 1 ^)) dl ■ ■ ■ {y-a^x^Y-iy-a^x^r • • • (y-a^x 1 ^)) 
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is the complete factorization of f(x,y) in C((x))*[y], with a^t) = J2 k Pi k t k , a 
series with real coeficients pi k , and etj(t) = Cj k t k , series having at least one 
nonreal coefficient, then fi(x,y) = (y — ai(x 1 / <H )) di modulo x n+1 , y a, = r«/r. 
The code for factoriza is 

> factoriza: =proc(f,n) 

> local FF, RR, i, redu; 

> FF:=[f];RR:=[l,[l]]; 

> while (fiscal(FF)<>0) do 

> for i from 1 to nops(FF) do 

> if reduccion(FF[i],n)[2][l]<>l then 

> redu:=reduccion(FF[i],n);RR:=cambio(RR,i,redu[l],nops(redu[2])); 

> if i < nops(FF) then FF:=[op(FF[l..i-l]),op(redu[2]),op(FF[i+l..nops(FF)])]; 

> else FF:=[op(FF[l..i-l]),op(redu[2])]; end if; 

> else end if;end do; 

> end do; 

> RETURN( [factor (FF) ,RR] ) ; 

> end proc: 

The routine rotacion takes a polynomial f(x,y) and looks for an integer 
n > such that the substitution x = x + ny, y = —nx + y makes it quasi- 
monic. Afterwards it divides it by a nonzero constant making it monic. The 
algorithm produces a [g/c,n], where g is obtained from / via the substitution, 
c is the coefficient of the highest power of y in the polynomial g, and n is the 
integer that makes this substitution work. If f(x,y) is monic, the algorithm 
takes n = and hence it returns f(x,y) again. The code for this routine is 

> rotacion: = proc(f) 

> local c,d,g,n; for n from to infinity do 

> subs({x=x+n*y, y=-n*x+y},f); 

> g:=collect(simplify(%),y); d:=degree(g,y); 

> c:=coeff(g,y,d); 

> if degree(c,x)=0 then RETURN ( [g/c , n] ) ; end if; 

> end do; end proc: 

The routine revisor takes as entry a polynomial f(x,y) that is monic in 
y. If f(x,y) = g(x,y) n , then revisor return a polynomial g(x,y), if g(x,y) 
is linear in y, i.e. of the form y — a(x). Otherwise, it returns the empty set. 
The code for this routine is 

> revisor: =proc(f) 



19 



> local 1; 

> 1:= factors(f)[2]; 

> if degree(l[l][l],y)=l then RETURN(monico(l[l][l])); end if; 

> RETURN (); end proc: 

limite is the main routine and it is built out of the previous ones. It 
receives as entry polynomials f(x,y),g(x,y) (not necessarily monic) and an 
integer n > . It computes the quotient q — f/g- Then it computes 

h = y(g(df/dx) - f(dg/dx)) - x(g(df/dy) - f(dg/dy)). 

Then it performs an appropriate rotation hi = rotacion(h), so that hi is 
monic in y, and factors hi into irreducible polynomials hi = bl • ■ ■ bs. Next 
if applies factoriza, to an approximation of n > 0, and returns for each factor 
bi, a list L formed by two lists 

LI = [\pi(x,y),...,p n (x,y)] y L2 = [r, [ri,...,r n ]] 

such that r — r± • • • r n y bi(x r , y) = pi(x ri ,y) ■ ■ -p n (x rn ,y). 

After this, it computes the list S = [pi(x Tl , y), ...,p n (x Tn , y)]. Eachpj(x r % y) = 
(y — a i (x ri )) di ) (which would correspond to a real Puiseux series) or having 
the form pi = (u(x,y)) d % where u(x,y) is monic in y, and of degree larger 
than one (which would correspond to a nonreal Puiseux series). Observe that 
bi(x,y) = pi(x Tl/T ,y) ■ ■ ■ p n (x Tn/r ,y). 

After this, revisor is applied to S and returns a list 

H = [qn(x,y), q i2 (x,y), q ik (x,y)\, 

where th qij(x,y) — (y — ciij(x rij )) are exactly those elements of S , without 
the power g^, that correspond to real Puiseux series. That is 

qn(x,y) dil = p a (x ril , y), q ik (x, y) d * k =p ik (x nk ,y). 

In this way qn{x,y) = y- a a (x rn ), q ik (x, y) = y - a ik {x nk ). 

Then it constructs f\ and g±, the polynomial obtained by applying the 
substitution x = x + ny, y = —nx + y to / and g (using the same n found 
for the rotation of h). After this it makes the list T = [an(x ril ), a ik (x rtk )]. 
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Then it computes R = [aii(x ni ^ r ), a,ik(x rik / r )], and then it computes the 
list P of pairs 

P = [[A(x,a a (:r r<l/r )),2i^ 

excluding those trajectories not passing through the origin ( aik(0 Tii ^ r ) ^ ) 
and therefore the list of limits 

Q\i] = [limj/^x, a^^/^/ftlx, au^'^))), .., limj/^x, a ife ^/ r ))Ai(^ %(^ /r )))L 
for i — 1, s. 

The limit exists if all the values of the list Q[i], for all bi, % — 1, s, are 
equal. 

> limite:=proc(f,g,n) 

> local F,G,a,b,Q,e,P,k,T,j,N,fl,gl,i,S,H,L,hl,h,q,R,t; 

> q:=f/g; 

> h: = simplify (expand(y* (g*diff(f ,x)-f*diff(g,x) )-x* (g*diff(f ,y)-f*diff(g,y) ) ) ) ; 
if h=0 then RETURN (q) end if; 

> hl:=rotacion(h); 

>F:=factors(hl[l])[2]; G:=[]; for b from 1 to nops(F) do G: = [op(G),monico(F[b][l])]; 
end do; 

> for a from 1 to nops(G) do Q[a]: = []; 

> L:=convert(evalf(factoriza(G [a], n)), rational); S: = []; t : =root(x,L [2] [1], symbolic); 

> for i from 1 to nops(L[l]) do 

> S:=expand([op(S),subs(x=x"L[2][2][i],L[l][i])]); end do; 

> H:=map(revisor,S); N:=hl[2]; 

> subs({x=x+N*y, y=-N*x+y},f);fl:=collect(%,y);subs({x=x+N*y, y=- 
N*x+y},g); 

> gl:=coUect(%,y); T:=Q; 

> for j from 1 to nops(H) do 

> T:=[op(T) r coeff(H[j],y,0)/coeff(H[j],y,l)]; end do; R:=subs(x=t,T);P:=Q; 

> for k from 1 to nops(R) do if subs(x=0,R[k])=0 then 

> P:=[op(P),subs({y=R[k]},[fl,gl])]; end if; end do; 

> for e from 1 to nops(P) do 

> Q[a]:=[op(Q[a]),Umit(evaJf( P[e][l]/P[e][2] ),x=0)]; end do; 

> end do; [seq(op(Q[a]),a=l..nops(F))]; 

> end proc: 
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4 Calculation 

Next, we give some examples illustrating the performance of the algorithm: 

1. limite(6*x"3*y,2*x"4+y"4,20); [2.033104508, -2.033104508, 0.]. Con- 
sequently, lim( Xj j / )_ > ( 0j o) 6x 3 y/(2x 4 + y 4 ) does not exist. 

2. limite((x"3+y"3),(x"2+x*y+y"2),20); [0.,0.,0.]. In this case lim (;Ei ,, H(0i0) (a; 3 + 
y 3 )/(x 2 + xy + y 2 ) exists and equals 0. 

3. limite(6*x"3*y,2*x"4+y"4,20); [2.033104508, -2.033104508, 0.]. The 
limit lim( Xj j / )_ > (o j o) 6x 3 y/(2x 4 + y A ) does not exist. 

In the following the limit exists only in 4. and 5. 

4. limite((x"4-y"2+3*x"2*y-x"2),(x"2+y"2),20); [-1., -1., -1.]. 

5. limite( (x"2-y"2), (x"2+y"2),10); [1, -1.]. 

6. limite((x"6-y"4+3*x"2*y"3-x"4*y),(x"4+y"4+x"2+y"2),20); [0., 0., 0., 
0.] 

7. limite(x,x"2+y"2,30); [undefined] 

8. limite(y"4,x"4+3*y"4,50); [0.3333333333, 0.] 

9. limite(6*x"3*y,2*x"4+y"4,10); [0., 2.033104508, -2.033104508]. 

10. limite(x"4*y"4, (x"8+y"8)"3, 20); [0., 0., Float (infinity), Float (infinity)]. 

5 Discussion and Conclusions 

The algorithm developed in this article provides a method for computing 
limits of quotients of real polynomials in two variables which has proven 
to be more powerful in handling these type of limits than other existing 
algorithms. The following examples compare the performance of the routine 
here presented and that of Maple 12. 

1. Example number 5 in the previous section shows that lim^j^^o)^ 4 — 
y 2 + 3x 2 y— x 2 ) / (x 2 +y 2 ) is equal to —1. However, Maple 12 is uncapable 
to compute it. 
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2. Example 8 shows that lim( a . iJ/ )_ > ( ,o) 2/ 4 / (^ 4 + % 4 ) does not exist, an 
example that Maple 12 computes successfully. 

3. \im iXty) ^( 0fi) (x 6 -y A + 3x 2 y-x' i y)/(x' i + y 4 + x 2 + y 2 ) is equal to zero, 
a problem that defeats Maples ' routine. 

The theoretical method for computing limits developed in this article 
applies to quotients of two real analytic functions. However, the algorithm 
was only implemented for polynomials. The next logical step would be to 
extend this algorithm to cover this general case. 

In a sequel article we will develop a more general method dealing with 
limits of quotiens of real analytic functions in several variables. 
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